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ABSTRACT 


The fundamental methods are described for the general 
spacecraft trajectory design and optimization software 
system called Copernicus. The methods rely on a unified 
framework that is used to model, design, and optimize 
spacecraft trajectories that may operate in complex gravi- 
tational force fields, use multiple propulsion systems, and 
involve multiple spacecraft. The trajectory model, with 
its associated equations of motion and maneuver models, 
are discussed. 

Key words: Copernicus, Spacecraft Trajectory Optimiza- 
tion Software. 


1. INTRODUCTION 


Copernicus [17] is a general trajectory design and opti- 
mization system intended to solve a wide range of space- 
craft trajectory problems in a robust and efficient man- 
ner. It facilitates the design and optimization of simple 
to complex spacecraft trajectories. Some examples from 
recent Project Constellation studies include: trans-lunar 
[3], trans-Earth [16] and trans-lunar abort return trajec- 
tories [14]. Interplanetary mission applications include: 
design of a 10 year 32-asteroid tour using a low-thrust 
propulsion system, Mars sample return missions and tra- 
jectory design for the Lunar Crater Observation and Sens- 
ing Satellite (LCROSS). 

The three main problems addressed by the system are: 

1. Modelling and Open Loop Simulation: The system 
can simulate open loop trajectories in any force field, 
using one or more spacecraft, and including impul- 
sive and/or finite bum maneuvers. 

2. General Targeting via Nonlinear Root Finding: The 
system can search for the values of the independent 

variables required to satisfy a set of constraint func- 
tions. This is required for trajectory targeting prob- 
lems, such as orbital boundary value problems. But 


more generally, it is possible to formulate and solve 
general systems of nonlinear equations that depend 
on a preselected set of independent variables. 

3. Optimization: The system can extremize the value 
of a function consistent with the constraints. Any 
computable function can be extremized. Typically, 
the optimization mode of the system is used for ma- 
neuver design and optimization. This includes de- 
termining the times and parameters that describe the 
maneuvers, which can be impulsive or finite bums. 

A key component of the system is the trajectory model 
that describes how a complete trajectory for a single or 
set of spacecraft is modelled. The supporting compo- 
nents of the trajectory model include the equations of mo- 
tion, the propulsion system, the control parametrization, 
the independent variables, and the constraint functions. 
This paper describes and summarizes some of the theo- 
retical aspects of these components. Some of these along 
with the details of the system architecture have been doc- 
umented previously as the development of the system has 
evolved over the years[8, 9, 10]. A discussion on future 
extensions and proposed formulations for the finite bum 
maneuver modeling is also given. 


2. TRAJECTORY MODEL 

The trajectory model uses a building block called a seg- 
ment. A segment is a trajectory arc that can have veloc- 
ity impulses and/or a finite bum maneuver. Depending 
on how the relationship between different segments is de- 
fined, single or multiple spacecraft trajectories that may 
interact can be modeled. The generality of how segments 
are related to one another facilitates the modeling of com- 
plex missions. The simplest mission may have only one 
segment. Complex missions may have tens or hundreds 
of segments. 

The segment arc is bracketed by two node points. The 
node points are tagged with epochs that are referenced to 
a reference epoch, t epoc h . The epoch of the initial node 
is t 0 and the epoch of the final node is tf. There is no 



restriction on the values for to and tf (to = tf, to < tf, 
to > tf ). The segment duration is A t = tf -t Q . The 
segment nodes are uniquely defined by specifying any of 
the pairs (t 0 ,tf) , (to, At) , (A t,tf). 

The segment state vector x is 

x = ( r T v T m ) T . (1) 

The initial position vector is r 0 , the initial velocity vector 
is Vq , and the initial mass is m 0 “. Both v and m can 
have discontinuities at the initial node. For velocity, the 
discontinuity is an impulsive maneuver Avo* For mass, 
the discontinuities are due to the propellent mass con- 
sumed by the impulsive maneuver, and independent mass 
discontinuities before and after the impulsive maneuver. 
The independent mass discontinuities can be used to rep- 
resent acquisition or removal of other spacecraft, and/or 
stage masses. After these discontinuities are accounted 
for, the segment is propagated from to to tf if to ^ tf. 
At the tf node, velocity and mass have the same type of 
discontinuities as at the to node. The propagation direc- 
tion is considered as forward if t o < tf and backward 
if to > tf. If to = tf, then the propagation direction 
is specified to properly account for how the velocity and 
mass discontinuities are used. Let id = +1 for a forward 
propagation and id = - 1 for a backward propagation. 

Across either the t 0 or the tf node, position remains con- 
stant. The velocity at the end of the node is 

v + = v“ + (id) Av (2) 

The mass across the node evolves as 

m — > m _+ — > m + ~ — » m ++ 

Any one of these is specified, and the remaining three are 
determined by knowing the propagation direction, the in- 
dependent mass discontinuity Am" before the impulse, 
the independent mass discontinuity A m + after the im- 
pulse, the magnitude of the impulsive maneuver, and the 
exhaust velocity c of the engine doing the impulsive ma- 
neuver. Define the factor 

/ = e~ Av / c (3) 

There are two main cases: 

1. m or m~ + is specified; depending on which one 
is specified, the other is solved by using 


m + = m + Am 

(4) 

m +- = (m _+ ) (f) ld 

(5) 

m ++ = m H + Am + 

(6) 


2. ra + or m ++ is specified; depending on which one 
is specified, the other is solved by using 

m ++ = m + ~ + A ra + (7) 


then 

m~+ = (m + ~) ( f)~ %i (8) 

m = m + ~ — Am - (9) 

As stated before, the segment initial state vector is 

x o = ( *o" v 0 - T m 0 ““ ) T . (10) 

After accounting for the potential velocity and mass dis- 
continuities at the to node, the initial condition for the 
state vector that is propagated from to to tf (if to ± tf) 
is 

x 0=( r ([ V(J- T m^+ ) T . (11) 

The state vector after the propagation to tf is 

Xf = ( rj vJ T m,Q- ) . (12) 

The final value of the state vector after all possible veloc- 
ity and mass discontinuities are accounted for at the tf 
node is 

x+ = ( rj v+ T mr ) T . (13) 

Sets of segments are used to construct complete trajecto- 
ries for one or more spacecraft. Segments can be con- 
nected to other segments by forcing the complete state or 
a subset of it to inherit the complete state or a subset of it 
from either node of another segment. Segments can also 
be completely disconnected to model independent trajec- 
tories representing other spacecraft, for example. These 
segments can later be constrained to be connected via 
continuity constraints to other segments if needed. A 
mission consists of the set of all segments. 

A mission has a base frame which is required to be non- 
rotating and centered at a celestial body. All state trans- 
formations use the base frame as a hub through which all 
state, maneuver, and state function transformations are 
made. Each segment has a state input frame, a propa- 
gation frame, and a function output frame. The propa- 
gation frame is the frame used for the integration of the 
equations of motion. It is required to be non-rotating 
and centered at a celestial body. The initial position 
and velocity state data (r 0 , v^") is given in the state in- 
put frame. The function data for a segment includes all 
functions that can be computed from the distinct segment 
states (r 0 , r / , v^“ , Vq~ , , and so on) and are com- 

puted in the function output frame. The state input and 
function output frames are required to be centered at a 
celestial body and can be either fixed or rotating. Each 
segment is allowed to have its own state input, propaga- 
tion, and function output frame. A simple mission may 
have all of these frames be the same. 

The individual segment impulsive maneuvers are each 
referenced to an impulsive maneuver frame. The impul- 
sive maneuver frame can be either fixed (non-rotating) 
or osculating along the instantaneous trajectory path. In 



either case, the impulsive maneuver frame is referenced 
to a frame centered at a celestial body. A fixed impul- 
sive maneuver frame is typically the same as the seg- 
ment propagation frame. Alternatively, an example of 
a common osculating impulsive maneuver frame is one 
that uses the instantaneous velocity vector as a basis vec- 
tor. The remaining two basis vectors can be constructed, 
for example, by choosing the instantaneous angular mo- 
mentum vector, if it exists, as one of the basis vectors. 

A segment vector s l is defined for each segment i to be 
a mixed variable type vector containing at least all of the 
independent and dependent variables and functions asso- 
ciated with the segment. Additionally it contains all the 
frame definitions, how and what data is inherited from 
other segments, what functions need to be evaluated from 
the state vector at the distinct time tags, which variables 
are search variables, and so on. This list of information 
is large. As an example, assume that a segment is purely 
ballistic with impulsive maneuvers at both nodes. A sub- 
set of the segment vector s l is 

s i = f *o,rJ,Vo T ,mo - ,Amo , AvJ, AmJ, \ T 

\ tf,AmJ,AvJ,Amj,r f ,vj + ,m^ + , J 

(14) 

The complete mission with n segments can be computed 
uniquely in open loop (simulated) by processing all of the 
information contained in s i ( i = 1 

Every segment i has a variable vector s l x and a function 
vector Sj. s* is a subset of s l and contains only the vari- 
ables that have been identified to be independent search 
variables to be determined, is a subset of s l and con- 
tains only the functions that have been identified to be 
constrained or added to an objective function. 

The state vector for any segment is governed by the first- 
order vector equation of motion 


dt 


( 


r 

v 

m 


g( r , v,m,t, p) + ^u(i) 

mi 

c(t) 


(15) 


where g is the ballistic acceleration per unit mass, p is 
a vector of non-state or non-control parameters such as 
force model and vehicle parameters, T is the finite en- 
gine thrust force, u is the unit thrust direction, and c is 
the engine exhaust velocity. T, c, and u are possible con- 
trol functions. The information needed to evaluate these 
functions are also contained in s\ These equations are 
solved as an initial value problem from to and tf. The 
only allowable discontinuities along the integrated solu- 
tion are associated with the control terms and only if an 
optimal control formulation is used to determine the val- 
ues of the control functions. 


The natural ballistic acceleration vector g is expressed 
in the segment propagation frame which is centered on a 
main celestial body q, and is non rotating. An example 


expression for g is 


g(r,i) = - 


Gm c 


rib 


- g E 


m 


3 = i 


M \r~rM 3 


+ 


d(*)/ 


~^~^drag "t" &srp H" ^non— spherical 


where G is the universal constant of gravitation, m Cb is 
the mass of the main celestial body, rib is the number of 
additional celestial bodies, and rrij is the mass of celestial 
body j . It is assumed that the positions of the rib celestial 
bodies is known as an explicit function of time with re- 
spect to Cb. Depending on the problem, additional terms 
accounting for other common accelerations, such as at- 
mospheric drag(ad r a ff ), solar radiation pressure(a srp ), 
and non-spherical celestial bodies(a non spherical), are 
added to the vector function g. The equations of motion 
are integrated numerically with an explicit numerical in- 
tegration algorithm. The trajectory arc that connects the 
to node to the tf node is either ballistic or it is a finite 
bum arc. If it is a finite bum arc, one of two main finite 
bum models are used. These are described below. 


3. FINITE BURN MANEUVER MODELS 


The modelling of a finite bum arc is based on one of two 
possible models. The first model is the parameter model 
that uses a finite set of parameters to describe the time 
evolution of the control variables. This finite set of pa- 
rameters form part of s^. The second model is the op- 
timal control model which requires augmenting the state 
vector to include the costate vector which is the vector 
of Lagrange multipliers adjoined to the state. The ini- 
tial conditions of the costate vector, or the initial values 
of a transformation vector used to uniquely determine the 
initial costate vector, form part of s* . 


3.1. Simple Engine Model 

The following discussion refers to the trajectory arc for a 
particular segment between to and tf (to < t < tf ) . The 
finite bum acceleration is given by die thrust acceleration 
term 


thrust ~ 

Tit ) , . 
: a\ U W 

m(t) 

(17) 

The mass rate is 

m = 

-m 

c(t) 

(18) 

and the minimal set of control constraints are 


T • < 

i min „ 

T(t) < T max 

(19) 

|u(t)| = 

1 . 

(20) 

c(i) > 

0 

(21) 



The type of engine being modelled determines which of 
these are the controls. A constant exhaust velocity en- 
gine that can only be turned on and off with T = T max 
when it is on has c(t) = c and T is piecewise constant 
(0 or T max ). If the engine can be throttled then T(t) 
needs to be determined consistent with its lower and up- 
per bounds. A simple model assumes that T(t) = T, a 
constant, with 0 < T < T max . 


An engine that is limited by the input power constrains 
T(t) and c(t) to satisfy 

P(t) = \nt)c{t) (22) 

with 

Pmin ^ m < ^max (23) 

along with bound constraints on either T(t) 


Pmin < T(t) < T max (24) 


or c(t) 

Cm in ^ c(f) ^ Cmax (25) 

In this model, any pair (T(t),P(t)), (T(t),c(t)), 
(c(t) , P(f )) are the control variables. The power P can 
be obtained by an internal source or by external source 
such as a solar electric system and will depend in part on 
the radial distance from the Sun. Optimal control theory 
has been used to determine the values of the controls that 
affect directly the thrust acceleration magnitude, namely 
T, c, P, for this system[9]. 


3.2. Detailed Engine Model 


For the detailed engine model the controls can be Thrust 
(T), specific impulse ( I sp ), exhaust velocity (c). Power 
(P), mass flow rate (m), Thrust acceleration (T acc ), or 
Thrust over weight ( T/W ). In addition, the efficiency 
( 77 ) of the engine can be modeled which relates the in- 
put power ( Pinput ) to the power available to the engine. 
Any valid combination of the engine controls can be se- 
lected (e.g., Thrust and 7 sp , Power and c, etc.) Any of 
the control parameters can be specified as constants, or 
can be modeled using selected engine control laws. The 
selected engine controls must satisfy the equations: 


where BB and DD are specified coefficients. Alternately, 
polynomial expressions for mass flow rate and thrust as a 
function of power can be used: 

m = Cf Ml) + c m (2) • P + M 2 ) • -P 2 (33) 
+c m (4)-P 3 +c m (5)-P 4 ] 

T = c r [c t (l)+c t (2)-P + ct(2)-P 2 (34) 

+c t (4)-P 3 +c t (5)-P 4 ] 

Where Cf and c t are coefficient arrays for modeling LILT 
effects and solar panel degradation, respectively; Cf is the 
duty cycle fraction. 


The engine can use an internal power source, or be mod- 
eled as a solar electric engine. For a solar electric engine, 
a custom solar panel power model can be used: 


1 + g(4) ■ r + g{ 5) • r 2 r 2 
•[t(l) + f( 2 ) • e^ 3 ^' At + f(4) • Af] 


Where g and t are coefficient arrays, At is the accumu- 
lated time the vehicle has been in space, Po is the refer- 
ence power at 1 AU at the start of the mission, r is the 
distance from the Sun, and 0 is the fraction of the Sun 
that is visible. A variety of eclipse models can be used, 
from a simple cylindrical model to a complex model that 
includes the umbra, penumbra, and antumbra. The ac- 
cumulated time in space can be tracked across multiple 
segments, allowing for multiple coast and thrust arcs in 
the mission. 


The remaining discussion is limited to discussing the 
computation of the thrust direction unit vector u (t). 


3.3. Parameter Steering Model 


The thrust direction unit vector is parameterized by two 
spherical angles referred to the finite bum control frame. 
This frame can be the same as the segment propagation 
frame, but it can also be an osculating frame such as one 
that uses the instantaneous velocity vector as a basis di- 
rection. In terms of the spherical angles, 


V = 

P 

Pinput 

(26) 

p = 

i. T . c 

(27) 

fri = 

-T/c 

(28) 

c = 

isp ‘ go 

(29) 

Tacc = 

T/m 

(30) 

T/W = 

T/{m ■ g 0 ) 

(31) 


The efficiency of an electric propulsion systems can be 
modeled as a function of specific impulse using the equa- 
tion: 


_ BBxP sp 
V I 2 p + DD 2 


( 32 ) 


( cosa(i)cos/3(i) \ 

sin a (t) cos /3 (t) | (36) 

sin /3(f) ) 

The angles are time functions of the form 

a(t) = a 0 + d 0 (t - t 0 ) + ^ ^ 

+a a sin (cj a (t - t Q ) + <f> a ) (37) 

Pit) = A) + Po(t ~ *o) + ^ 2 

+ap sin (cjp (t - to) + 4>p) (38) 

where (to <t <tf). Each of the constants in Eqs. 37- 
38 can be part of s* . As defined, the functions for the 



For the force model considered, 


spherical angles a and 0 admit constant, linear, quadratic, 
and sinusoidal terms. For most practical applications, the 
constant and linear terms are sufficient. For longer fi- 
nite bum durations, both the quadratic and the sinusoidal 
terms may be needed. 

A more general model that attempts to emulate the po- 
tential complexity that the steering of a finite bum can 
exhibit is referred to as the Single Axis Rotation (SAR) 
Model. In this model a reference unit vector u re f is ro- 
tated about a rotation axis n through a time varying angle 
7 . The reference unit vector u re f is referenced to a con- 
trol frame that in turn is defined by basis vectors defined 
in a frame centered at a given celestial body c&. The 
spherical angles that define the rotation vector n and the 
rotation angle 7 are functions of time similar to the time 
functions given by Eqs. 37-38. A special case exists if 
u (t) is confined to a plane normal to the rotation axis n; 
here the new reference unit vector u re j that gets rotated 
is 

/ _ (n x u re /) x n 

Uref |(nxu re /)xn| 

Depending on the parameters used to define u re /, n and 
7 , u(t) can be constant, or it can rotate at a constant rate 
in plane normal to n, or it can precess and nutate relative 
to a cone whose central axis is the rotation axis n. In 
general, the unit vector thrust direction is calculated: 


f(x,t,u c )= g(r, v,m, i, p) + ~~u(t) (44) 

\ j 

\ c(t ) / 

the Hamiltonian is 


H = ( g(r, v,m, t , p) + — u J +A„ 


and the costate vector equations are 


dr) Av 


(£) 




dm v V dm m 2 * 1 


Extremization of H with respect to the thrust unit vector 
results in the well known thrust unit vector steering law[5, 

<49> 

and is referred to as the Primer Vector. 


U (t) = U re f(t) T n(t) (1 - COS Tl(t) 

+ C 0 S 7 (t)u re /(f) +sin 7 (t) (n (t) x u re f(t)) 


The solution requires knowing the costate vector A(to); 
and it forms part of s^,. An estimate for A (to) is ob- 
tained by using an adjoint-control transformation^]. If 
an estimate for uq is available, 




^ Vo A<y 0 U.O 

(50) 

3.4. Optimal Control Steering Model 

where uq 

/ COS Qio COS 0 q \ 




Uo = I sin a 0 cos 0 Q I 

(51) 

In the Mayer form of the optimal control problem[l, 4] 


V sin 0 q J 



where the cost function to maximize is of the form 

J = (39) 

the optimal control Hamiltonian is 

H = A T f (40) 

where A is the costate vector of Lagrange multipliers ad- 
joined to the state vector and f are the state equations, 

x = f(x,t,u c ) (41) 

and u c is a vector of control variables. A first-order nec- 
essary condition is the vector differential equation for A 


Additionally, the control vector u c is chosen to extremize 
H at all points on a solution consistent with the control 
constraints if they are present. The explicit costate vector 
is 


The time rate of change of \ Vo is 

= “I" A Vo Uo 

where 

( —d 0 sin a 0 cos 0 Q — 0 O cos a 0 sin 0 Q 
do cos Q!o cos 0 Q — 0 Q sinao sin/3 0 

ko COS P 0 

Using Eq. 47 


A-Uq + A-y 0 I 1 U.Q A-y 0 U.Q ■ (54) 


The basic transformation is summarized as 
/ A<y n \ / A ?Jn \ 


Au 0 
A^o 
<X(hPo 
& <hPo 
A rn.n 


■ ( 55 ) 



The necessity of explicitly specifying all of the compo- 
nents of the costate vector has been removed in favor of 
specifying the four control related variables (a, /3, a, /3) 
and three costate related variables (A v , \ v , A m ). Without 
any prior knowledge of the solution for u(t), obtaining 
an estimate for these remaining quantities is possible, but 
is problem dependent and requires further knowledge as- 
sociated with the additional transversality conditions as- 
sociated with the optimal control problem [9]. If a pa- 
rameter model solution exists, a method to obtain these 
for a constant exhaust velocity engine is discussed in the 
last section of this paper. The use of the transformation 
replaces the explicit costate vector in s^, with the adjoint 
control variables. 

It is important to note that even though an optimal control 
formulation leads to a well defined multi-point boundary 
value problem with functions that include transversality 
conditions on the states and costates at distinct times, 
these conditions are not used in the system. This is in- 
tentional. Instead, the method used is best described as 
a hybrid method where some of the unknown variables, 
which include the initial values of the costates or func- 
tions of them, form part of the variable vector and the 
cost function is optimized directly with an optimization 
algorithm. The advantage being that the control follows 
the control optimality condition, without explicitly deriv- 
ing and enforcing the transversality conditions. 


4. TARGETING AND OPTIMIZATION 


Let x p be the n p x 1 vector that contains the set of vari- 
ables to be determined. x p is an augmented vector that 
contains all of the segment variable vectors 4(* = 1 - n )- 
The function vector is the n c x 1 function vector c that 
contains all of the segment function vectors = 

c is functionally dependent on x p , 

c = c(xp) (56) 

The constraints on c are 

c l < c(xp) < c u (57) 

where c l and c u are constant vectors of lower and up- 
per bounds, respectively. An equality constraint on an 
element j of c is defined by setting 

<4= C " = S* (58) 

where c* is the required target value. Simple one sided 
inequality constraints Cj, are set by choosing 

— oo < Cj < c 'j or c l j < Cj < -boo (59) 

If rip = n c and all of the constraints on c are equality 
constraints, then the problem requires solving a nonlinear 
system of equations. If n p < n c and all of the constraints 
on c are equality constraints, then a minimax solution is 


needed where it is necessary to find the local minimum 
of the function 


max 


[c(x P )L 


, j = 1 ...n c . 


(60) 


If n p > n c , the system is underdetermined, the minimum 
of the function in Eq. 60 zero, and a solution is the min- 
imum norm solution from the initial value of x p . For an 
optimization problem, let the element k of c, c*., be the 
objective function with the condition that c l k ^ cj£. The 
nonlinear constrained optimization problem is to mini- 
mize or maximize Cfc subject to the constraints given by 
Eq. 57. The problem has been cast in a standard form 
currently used by state of the art nonlinear constrained 
optimization codes such as SNOPT[7, 15]. 

Regardless of whether a targeting or optimization prob- 
lem is required to be solved, it is desirable to at least 
estimate as accurately as possible the n p x n c Jacobian 
matrix dc/dx p . Numerical finite difference methods are 
the most common methods for this purpose, even though 
they have inherit difficulties especially for the functions 
in c that are highly sensitive to the elements in x p . For 
example, the (i,j) element of <9c/dx p can be approxi- 
mated using a central difference approximation 

~ c i{ x Pj ~ c i{ x Pj — ^ x Pj) ( 61 ) 

dx Pj ~ 2Ax Pj 


where Ax Pj is the positive perturbation stepsize for the 
x Pj element of x p . A very accurate method to compute 
<9c/<9x p uses state transition matrices and analytical gra- 
dient expressions for all of the functions and the states 
required required to compute those functions[18, 11, 12]; 
however, this method is not practical for a generalized 
system because of the overhead required to derive all of 
the required relationships between the large number of 
search variables and constraint functions. 

Examples that describe how to cast mission design and 
trajectory optimization problems into a targeting or op- 
timization solution method for both impulsive and finite 
bum maneuvers are described in[13]. 


5. FUTURE CONSIDERATIONS 


One of the areas currently being examined is an alter- 
nate formulation for the finite bum thrust steering model. 
For the parameter basic steering model, assuming only at 
most a quadratic function for the spherical angles used to 
represent u(t), an augmented state vector that includes 
the parameters that define the thrust unit vector u(t) and 
its rate u(t) is being considered, 


x = ( r T v T 


P a p ) T (62) 


m a 



with 


Using S' (t 0 ) = 0, 


x(t,x) 


/ v \ 

g(r,v,m,t,p) + ^u(t) 
T(t) 
c(t) 

a 


p 

a 

p J 


(63) 


In this model, the initial conditions for the spherical an- 
gles and their rates form part of the segment variable vec- 
tor. With this model, a procedure is being investigated to 
automate the conversion of a parameter model finite bum 
solution to a finite bum solution that uses the primer vec- 
tor for the thrust unit vector. As an example, assuming 
that g only depends on r and t, it can be shown that an 
augmented state vector of the form 


<7i> 

The time history of X v (t) is assumed to be quadratic 
with negative curvature (A v < 0) and with a maximum 
value near the middle of the bum. If t m is the time of 
the midpoint of the bum then as an approximation set 
X v ( tm ) = 0. Having a parameter model finite bum solu- 
tion, Eq. 67 is used at t m to get an approximation for \ v . 
X v is assumed to be constant; this along with the condi- 
tions X v (to) = 1, X v (tm) = 0 the expression for A v (to) 
is 

K{to) = “A v (tm) [(u T G) (u) + U T u] t ( t m - t 0 ) 

(72) 

where X v (t m ) is the value of X v at t m and is given by, 


x = ( r T v T m a (3 a f3 X v X v A m ) T 

(64) 

has the following first order equation of motion 


x(i,x) 


V \ 

g( r , p) + £u 

_ T 
c 

a 

p 

a 


"p 

Xy 

A y 

V m?^ v J 


(65) 


where a , p are determined uniquely from u and u; and 
a, p are determined uniquely from u, li and ii; and 


u = Gu-{[(u T G)(u)] +u T u}u-2^u (66) 

Xy 


X y = Xy [(u T G) (u) + U T u] 

where 

G= dg (M) 

dr 


(67) 

( 68 ) 


With this formulation, an automated procedure is being 
developed that transforms a parameter model finite bum 
solution to an approximation of the nearby optimal con- 
trol finite bum solution. For a constant exhaust velocity 
engine use is made of the associated switching function 
S that determines the value for T, 



m c 


(69) 


The cost function is minimization of the propellant con- 
sumed by the finite bum. With this, S (to) = 0 and to is 
forced to be the switching time from zero thrust to T max . 
The scale for the costates is set by choosing 


X v (to) = 1 (70) 


1 

“ 1 + I [(u T G) (u) + u T u] tm (t m -t 0 ) 2 

(73) 

This procedure provides the initial estimate for all of the 
adjoint control variables needed in the augmented state 
vector, Eq. 64, at to . 

Another area being considered for possible inclusion in 
Copernicus is a general architecture for guidance algo- 
rithms and control laws. Current finite bum models 
in Copernicus can be complemented with guidance al- 
gorithms for specific applications. Finding a general 
method for the implementation of such algorithms is still 
under research. The implementation of a general relative 
motion capability between segments is also under con- 
sideration. This will allow a more flexible modeling and 
optimization of spacecraft rendezvous and formation fly- 
ing problems. Finally, a more complete vehicle model 
within the segment that incorporates the ability to model 
a spacecraft as a system of multiple masses, e.g.: tanks, 
consumables, dry mass, stages, propellant, etc. is being 
designed. 
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